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"*^ , ABSTRACT 

(N 

^S] ' We present a study of the circumstellar environment of IRAS 04158-1-2805 based on multi-wavelength observations 

and models. Images in the optical and near-infrared, a polarisation map in the optical, and mid-infrared spectra were 
obtained with VLT-FORSl, CFHT-IR, and Spitzer-IRS. Additionally we used an X-ray spectrum observed with 

i-J^ ■ Chandra. We interpret the observations in terms of a central star surrounded by an axisymmetric circumstellar disc, 

t*H' but without an envelope, to test the validity of this simple geometry. We estimate the structural properties of the disc 

and its gas and dust content. We modelled the dust disc with a 3D continuum radiative transfer code, MCFOST, based 
on a Monte-Carlo method that provides synthetic scattered light images and polarisation maps, as well as spectral 
energy distributions. We find that the disc images and spectral energy distribution narrowly constrain many of the disc 
!!* ' model parameters, such as a total dust mass of 1.0 — 1.75 • 10~* Mq and an inclination of 62° — 63°. The maximum 

grain size required to fit all available data is of the order of 1.6 — 2.8 /xm although the upper end of this range is loosely 
constrained. The observed optical polarisation map is reproduced well by the same disc model, suggesting that the 
geometry we find is adequate and the optical properties are representative of the visible dust content. We compare the 
^ ' inferred dust column density to the gas column density derived from the X-ray spectrum and find a gas-to-dust ratio 

fT^ ' along the line of sight that is consistent with the ISM value. To our knowledge, this measurement is the first to directly 

00 . compare dust and gas column densities in a protoplanetary disc. 
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1. Introduction However, the properties of discs around more massive 

stars and, of concern here, around the lower mass late M 

Accretion discs are key elements in star and planet for- dwarfs and brown dwarfs remain poorly known, because 

mation scenarios. They provide the material for accretion images for these discs are still extremely rare. As a conse- 

leading to star and planet building, they provide the energy quence, our knowledge of the circumstellar environment of 

^ . and material for launching jets, and they are the medium these objects is based solely on spectral energy distribution 

through which angular momentum is transported away. (SED) fitting. 

Knowing their geometrical and physical properties is im- ^^ ^^^ ^^ ^^^^^ ^ ^^ ^^^^ 04158+2805, 

portant tor understanding these processes and their evolu- , , , i , , ,, , , mi 

i;. o J- g^ low-mass star located near the substellar boundary, ihe 
Lion 

classification of IRAS 04158+2805 varies in the literature, 

Large-scale surveys have been performed to search for ^^^h most authors agreeing on a late spectral type, implying 

discs around young low-mass pre-main sequence stars (e.g., ^ very low stellar mass. The recent paper by Luhman (2006) 

Stapelfeldt et al. 2003, Schneider et al. 2005), the so-called concludes that its spectral type is consistent with a low 

T Tauri stars. A few tens of discs have been imaged and, j^^ss star close to the brown dwarf limit. IRAS 04158+2805 

for many of them, images are available over a progressively jg located at a distance of 140pc, in the L1495 east dark 

broader wavelength range, enabling deeper studies of the ^loud, which is part of the Taurus molecular cloud complex, 

disc properties. jt is surrounded by an extended reflection nebulosity seen in 

scattered light. It propels a well-coUimated, ionised atomic 

~ ; ~ ~~ 7 J A Aj r^i -1 iet seen in Hq, that extends at least out to 60 arcsec to 

bend ojfpnnt requests to: A. M. Giauser, e-maii: ,, ,, mi i • , ■ i i i i , i • ,i r i 

glauserOphys ethz ch ^^^ north. The object is probably located m the foreground 

* Based on "data' collected at the Canada-France-Hawaii of the large reflection nebulosity illuminated by V892 Tau, 

Telescope. The CFHT corporation is funded by the governments one of the rare Herbig Ae stars of the Taurus cloud, because 

of Canada and France, and by the University of Hawaii. Based the disc of IRAS 04158+2805 appears in silhouette over the 

also on data collected at ESO/VLT during observation program reflection nebulosity. This is the only such case in Taurus 

68-C.0171. to our knowledge. 
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In this paper, we apply an axisymnietric, inclined-disc 
model to our data, fitting near-infrared images, near- and 
mid-infrared photometry, mid-infrared spectroscopy, mil- 
limeter photometry, and near-infrared polarimetry to de- 
rive disc geometry, disc mass, and dust composition. The 
large size of the disc derived from the near-infrared images 
(radius of approximately 8 arcseconds in the I band, corre- 
sponding to 1120 AU at a distance of 140 pc) is intriguing, 
especially around a star of such a low mass. We will show 
that the data are compatible with an inclined disc without 
any need to invoke a larger envelope. The paper is organ- 
ised as follows. Section [2] presents the observations which 
are further discussed in section [31 In section [H we discuss 
the modelling of the dust disc. We estimate the structural 
model parameters of the circumstellar disc and we discuss 
the quality and the uniqueness of the solution found in sec- 
tion [S] We present our conclusions in section [Bj 



2. Observations and data reduction 

2.1. Optical imaging and polarimetry 

2.1.1. Observations 

IRAS 04158-^2805 was observed on December 12, 2001 
in the I-band with the FORSl/IPOL instrument. The 
weather conditions were good and the seeing was measured 
between 0'.'9 and 1" over the observation period. The total 
field of view (FOV) of FORSl/IPOL is 6f8 x 6.'8 in the 
Standard Resolution mode with a focal scale of O'.'2/pixel. 
Polarimetry was performed by inserting a WoUaston prism 
in the beam. The prism splits the incident light beam 
into two separate beams of orthogonal polarisation states, 
the so-called ordinary (o) and extraordinary (e) beams. A 
stepped half-wave plate retarder was placed at the entrance 
of the incident beam and was rotated by steps of 22.5 de- 
grees. The separation of the two o— and e— beams on the 
CCD is performed by the WoUaston prism and overlap of 
the two beams is avoided by inserting a 9-slit focal mask. 
Each slit in the mask provided a ^ 20" x 6.'8 FOV. Images 
of IRAS 04158-f 2805 were recorded at 8 retarder positions 
with an integration time of 3 minutes per frame. The im- 
ages were then combined to yield the Stokes parameters I, 
Q and U. 

2.1.2. Data reduction pipeline 

A dedicated data reduction pipeline was written using 
NOAO/IRAF (see, e.g., Monin et al. 2006). The images 
are first corrected for bias and bad pixels, and then flat- 
fielded. In the next step, the images went through a po- 
larisation extraction routine in which the normalised flux 
difference between the ordinary and extraordinary images 
is calculated for every pixel of the image, and a Fourier 
series was computed to derive the Stokes parameter I, Q 
and ifl: We first computed the normalised flux difference 
F from the target flux f: 



F{ 






(1) 



Thereafter we computed Q and U: 



N-l N-l 

Y^-F{e,) cos{Ae,) t/-^-n0,)sin(40,) (2) 



j=0 



i=0 



^ See the FORS user manual at |http://www.eso.org| The po- 
larisation level, P, is obtained by calculating P = ^JQ^ + V^ / 1 
and the position angle, O, by calculating O = l/2arctan((7/(3). 



The absolute errors were estimated by using two indepen- 
dent methods: flrst, from the photon noise on the e— and 
o— beams separately, and then propagating the errors in 
the calculations of Q, U, P and 6; second, by measuring 
the standard deviation on the 8 images from the half-wave 
plate rotation. Both methods gave consistent results, of or- 
der 0.3%. The final intensity map is presented in Fig.[TJ 



2.1.3. Instrumental polarisation 

We estimated carefully the instrumental polarisation at the 
centre of the FORSl field by measuring nearby (i.e., high 
proper motion) unpolarised targets. We observed GJ 781.1 
and GJ 2147, two high proper motion stars. Because the 
immediate solar neighbourhood is remarkably devoid of 
dust (e.g., Leroy 1993, 1999) the interstellar polarisation of 
nearby stars can be considered null. The average of 4 mea- 
surements on both GJ objects gives Pjnst = 0.02% ±0.03%. 
We therefore believe that FORSl/IPOL instrumental po- 
larisation is very low on-axis, well below 0.1% at the centre 
of the field, and we did not attempt to remove it from the 
measurements. 

On the other hand, Patat & Romaniello (2006) showed 
that FORSl presents a spatially variable instrumental po- 
larisation component. This component follows a radial pat- 
tern with an intensity scaling as 0.06 r^ (in % if r is in ar- 
cmin), ranging from 0.1% in the central region, one arcmin 
in radius, to 1% at the edge of the FOV. First, the absolute 
value of the instrumental polarisation is low (0.23% at the 
position of the source) compared to the observed polarisa- 
tion of IRAS 04158-^2805 which is always higher than 3%. 
Second, since the source is small compared to the FOV, we 
decided to neglect the spatial variation of the instrumen- 
tal polarisation on the FORSl detector because it does not 
vary significantly across the object. 



2.2. Near-infrared imaging 

On October 29 and 30, 2001, we used the near-infrared 
CFHT-IR camera (Starr et al. 2000) at the Canada-France- 
Hawaii Telescope to obtain H— and K—hand images of 
IRAS 04158-^2805 with a pixel scale of O'.'211/pixel and a 
total FOV of 3.'6. Conditions were non-photometric and the 
seeing during the observations was 0'.'65 at X— band and 
0'.'9 at iJ— band, as measured from the average FHWM of 
several unresolved point sources located in the FOV. With 
each filter, two series of 10 jittered images were obtained 
in separate sets. Each set of images was first reduced as an 
independent dataset in the following manner. All images 
were median-combined to create a sky frame, which was 
subtracted from each image. The images were then flat- 
fielded, registered based on the location of a bright point 
source in the field and median-combined. The two indepen- 
dent images per filter resulting from this procedure were 
then averaged together to produce the final images pre- 
sented in Fig. [TJ 
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Fig. 1. I—,H— and iiT— band contour plots of IRAS 04158+2805. The contour lines are on the levels /max • 2 " with 
n = 1, . . . , 9 or 10 for the iiT— band image respectively. 



2.3. Mid-infrared spectroscopy 

We use archival data of the Spitzer Infrared Spectrograph 
(IRS, see, e.g., Houck et al. 2004) observation from March 
4, 2004 (program request 3534848) which was done in the 
spectral mapping mode with the two low resolution chan- 
nels: Short-Low (SL; 5.2-14 ^m, A/AA - 90) and Long- 
Low (LL; 14-38 ^m, A/AA ~ 90). The mapping mode (3 
exposures across the target for each slit and 2 nodding po- 
sitions each) was chosen due to the small mispointing of 
Spitzer in the early mission. Therefore, the flux of the ob- 
ject is separated in two or three observations. To recover the 
photometric information, all three parallel exposures were 
summed for the SL while the overlap of the FOV of the 
LL-slit is such that the full flux is extractable by summing 
the first and last exposures. 

We used the final products from the Spitzer Science 
Center's IRS data-reduction pipeline (post-BCD). To al- 
low a background subtraction the two nodding observations 
were used to subtract them from each other. The data ex- 
traction was done with the Spitzer IRS Customer Extractor 
software (SPICE) provided by the Spitzer Science Center. 



2.4. X-ray spectroscopy 

IRAS 04158+2805 was serendipitously observed by the 
Chandra X-ray Observatory in a field pointing at the L1495 
East dark cloud around V410 Tauri. The observation was 
performed with ACIS-S on March 7, 2002 between 6:16 UT 
and 11:45 UT, with a useful total exposure time of approxi- 
mately 17700 s. IRAS 04158+2805 was located 11.'6 off-axis 
on the ACIS-S chip SI. This resulted in a rather distorted, 
extended point spread function (PSF) with considerable 
background contributions. The data were, however, suffi- 
ciently clean to extract a useful spectrum (see Fig. [2]). 

We reduced all ACIS-S data using standard analysis 
techniques following the Chandra CIAO Science Threadfl 
Specifically, we flagged bad pixels, applied CCD charge 
transfer inefficiency corrections to create a so-called level2 
events file. We extracted source counts from a circular 
area with a radius of 27'.'6 around the source. To define 
a background spectrum, a source- free circular area with 
a radius of 98'.'4 was extracted from the same chip. The 
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Fig. 2. X-Ray spectrum of IRAS04158+2805 observed with 
Chandra 



counts were binned into a source and a background spec- 
trum. Appropriate responses were created using the mkrmf 
task, and the ancillary file was obtained from mkarf. 



3. Results 

3.1. Imaging 

Contour plots of IRAS 04158+2805 at L, H- and K-bands 
are presented in Fig. [T] In the I-band, the object shows a 
bipolar refiection nebula geometry. The dark lane, tracing 
the equatorial plane, is seen in absorption over the back- 
ground light. It separates a prominent triangular nebulosity 
located to the North from a low surface brightness elon- 
gated counter nebula to the south. The counter nebula is 
undetected at H- and K-bands. At these wavelengths, the 
triangular shape of the main nebula is still visible, but its 
extension decreases with increasing wavelength. 

The maximum width of the northern nebula is mea- 
sured to be 15'.'8, 14'.'1 and 10'.'8 at 1-, H- and K-band, 
respectively. The northern nebula has a triangular shape 
whose opening angle is ^^130° at all wavelengths. Finally, 
the distance between the peak of the two nebulae at I-band 
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Fig. 3. Modeled SED (dashed hne) compared to observed 
fluxes: IRS spectrum (Sohd hne); Strom & Strom (1994, 
squares); Kenyon et al. (1990) and Kenyon & Hartmann 
(1995, circles); Luhman & Rieke (1998) and Luhman (2000, 
crosses); Motte & Andre (2001, diamonds); this paper, see 
sect. 13.21 f triangle) . The error bars of the photometric fluxes 
are smaller than their symbols and therefore not drawn. 



is 4.8±0.2 arcsec. In Sect. 14. 2[ we define additional mor- 
phological and photometric observables to compare models 
and observations. 



3.2. Aperture photometry 

To construct the spectral energy distribution of 
IRAS 04158+2805 shown in Fig. [3l we use mid-infrared 
IRAS flux measurements from Kenyon et al. (1990), optical 
and near-infrared photometry from Kenyon et al. (1990), 
Strom & Strom (1994) and from the 2MASS point source 
catalog. The 1.3mm continuum flux is from Motte & Andre 
(2001). Because IRAS 04158+2805 is extended, we suspect 
that the 2MASS photometry may be underestimated 
since these fluxes are listed in the 2MASS All-Sky Point 
Source Catalogue (PSC) and were obtained by PSF fitting. 
To check this point, we extracted the photometry from 
our own near-infrared images by using three point-like 
sources that have 2MASS photometry and use them as 
relative photometric standards. We obtained K — 11.04 
and H = 12.10 for IRAS 04158+2805, with estimated 
uncertainties of 0.05 mag. For the H-band this is 0.28 mag 
brighter than the 2MASS photometry and 0.13 mag for the 
K-band. While the brightness profile of IRAS 04158+2805 
is peaked in the near infrared and may account for this 
discrepancy, we note that from our images we find a 
difference of only 0.03 mag in the photometry when two 
apertures of 4'.'2 and 10'.'6 are used. Therefore, we attribute 
the difference in photometry of IRAS 04158+2804 between 
2MASS and our images to intrinsic variability of the source 
rather than an aperture size too small in 2MASS. 



3.3. Mid-infrared spectroscopy 

The Spitzer-IRS spectrum shows absorption features at 
~10/im and ^15.2/im. Both can be identified as silicates 



and carbon dioxyde ice, respectively (see Fig. [3]). Similar 
features are seen for example in low-mass protostars, where 
the cold envelope and/or outer disc cause the absorption 
(e.g., Watson et al. 2004). In IRAS 04158+2805, it is in- 
teresting to note that the silicate absorption is maximum 
slightly short of 10 microns, around 9.5 microns. This is an 
indication for small dust grains made of either amorphous 
or crystalline material such as enstatite (Schegerer et al. 
2006). We did not attempt to fit the exact shape of the 
feature and cannot conclude about the exact mineralogy of 
the grains. It is also interesting to find CO2 ice signatures 
in absorption in a source where the disc appears to domi- 
nate the SED (see section |4. 2. 3p . but the high inclination, 
nearly edge-on, probably provides the needed column den- 
sity along the line of sight. However, it is beyond the scope 
of this paper to discuss quantitatively the IRS spectrum. 
We will do so in a forthcoming paper. 

We extract the fluxes at 12/im and 25/im and compare 
with IRAS fluxes pubhshed in Kenyon et al. (1990). The 
difference is of the order of 0.5% at 12/im and 1.3% at 
25/im and is negligible therefore. 



3.4. I-Band polarimetric imaging 

The light from the northern photometric nebulosity centre 
has a linear polarisation of 3.3 ± 1.2%, averaged within a 
region of 2 arcsec (10 pixels) in diameter. Since neighbour- 
ing objects do not show a significant polarisation rate, this 
is an indication that the peak in the nebula is not arising 
from stellar photons seen directly but, instead, from scat- 
tered light. The polarisation is maximum at the corners 
of the triangular northern nebula and has values between 
35% and 41% (see Fig. 2]). For this determination we use the 
highest "reliable" value. We decided to use only those pix- 
els for which the polarisation rate topology shows a certain 
smooth behaviour. This was quantified by determining the 
standard deviation of the polarisation rate of each pixel and 
its 4 neighbours. If the standard deviation is higher than 
20% (typical noise level), the pixel is not used. 



3.5. X-rays 

The Chandra ACIS X-ray spectrum of IRAS 04158+2805 
is shown in Fig. [21 The spectrum was binned such that the 
minimum number of counts per bin was 10 before back- 
ground subtraction. The total number of counts in the 0.5-6 
keV range was 100, which allowed for only a rather simplis- 
tic spectral model that is, however, sufficient to derive use- 
ful estimates of the absorption column density to the X-ray 
source. The source spectrum is obviously very strongly ab- 
sorbed, with essentially no counts detected below 1.5 keV. 
The X-ray source is hard enough to produce counts up to 
6 keV in the observation. 

We fitted the spectral data in XSPEC (Arnaud et al. 
1996), using the vapec collisional-ionization equilibrium 
model and a photoelectric absorption component, the latter 
essentially being defined by the hydrogen column density 
Nh between observer and source. We first fitted the spec- 
trum with a single thermal component, assuming that all 
element abundances are at 0.3 times the solar values, refer- 
ring to the solar abundances of Anders & Grevesse (1989). 
This corresponds to the elemental composition typically 
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Fig. 4. Polarisation maps of the observed (left panel) and modelled {middle panel) nebula. The vector length for 100% 
polarisation is indicated by the bar in the lower right corner. Right panel: Comparison of the polarisation level as a 
function of position in the observed (solid) and modelled (dashed) nebula. These curves are estimated along the ridges 
(black curves) and symmetry axis (gray) of the nebula. 



measured in coronae of magnetically active stars (Giidel 
2004). 

We found a best-fit temperature of 5.8 keV or 6.7-10^ K, 
and a best-fit value for Nh of S.St^^ x 10^^ cm~^, where 
the error range refers to the 68% confidence limit. Given 
that the spectrum reveals only the hard end of the en- 
tire soft X-ray spectrum, the fit may be biased by fitting 
a hot component only. This bias may be relieved some- 
what by fitting a continuous distribution of emission mea- 
sures (EME0). We adopted an EMD with a prescribed 
shape such as found for less- absorbed T Tauri stars, de- 
scribed in more detail by Telleschi et al. (2007). The EMD 
model essentially consists of two power laws on each side 
of a peak. Neither the location of the EMD peak nor the 
power-law slope toward higher temperatures were well con- 
strained, but the results for Nh were robust and converged 
to Nh = S.eti.o X 10^^ cni-2. 

To summarise these two approaches, irrespective of 
the uncertainties in the intrinsic X-ray spectrum of 
IRAS 041584-2805, we find that the X-rays, which are pre- 
sumably emitted in the corona surrounding the central star, 
are attenuated by a gas column density of Nh ~ 3.5 x 
10^^ cm~^, with a one sigma range of (1.9— 6.1) x 10^^ cm~^. 



ronment. It includes multiple scattering, passive heating of 
the dust disc and thermal re-emission by the dust to pro- 
duce synthetic images in all four Stokes parameters at any 
wavelength, as well as SEDs. The dust temperature, the 
same for all grain sizes, is calculated assuming local ther- 
mal equilibrium. 



4.1.1. System geometry 



We consider a flared density structure with a Gaussian 
vertical profile p{r,z) — po{r) exp(— z^/2 /i(r)^), vahd for 
a well-mixed vertically isothermal, hydrostatic, non-self- 
gravitating disc. We use power-law distributions for the 
surface density S(r) = Sq (r/ro)" and for the scale height 
h{r) = ho (r/ro)^ where r is the radial coordinate in the 
equatorial plane and ho is the scale height at the radius 
ro = 50 AU. The disc extends from an inner cylindrical 
radius i?in to an outer limit i?out- The central star is rep- 
resented by a point-like, isotropic source of photons. 



4. Modelling of the dust disc 

In this section we explore the parameter space of a sim- 
ple model comprising the stellar photosphere and a power- 
law disc (see 14.1.11 below) to explain the observations of 
IRAS 04158-1-2805 and its associated reflection nebulosity. 



4.1. Modelling technique 

Synthetic images, polarisation maps and spectral energy 
distributions are computed using MCFOST, a 3D radiative 
transfer code based on the Monte-Carlo method. MCFOST 
is described in details in Pinte et al. (2006). MCFOST 
solves the full polarised radiative transfer in dusty envi- 



EMD defines the the emission measure EM—nennV as a 
function of the temperature, where ne,nH are electron and hy- 
drogen densities, respectively, and V is the volume. 



4.1.2. Dust properties 



We consider homogeneous spherical grains and we use the 
dielectric constants described by Mathis & Whiffen (1989) 
in their model A, which accounts for the interstellar ex- 
tinction law. Grain sizes are distributed according to a 
power-law n{a) ex a~^'^ with amin and amax being the 
minimum and maximum grain radii. The interstellar val- 
ues from Mathis & Whiffen (1989) are a,„in = 0.005/.tm 
and a,nax — 0.9/xm. The mean grain density is 0.5 g cm~'^ 
as a consequence of the high porosity (80%) of the grains. 
In this work, Omax is considered as a free model parameter 
to fit for. Extinction and scattering opacities, scattering 
phase functions and Mueller matrices are calculated using 
Mie theory. The dust and gas are assumed to be perfectly 
mixed and the grain properties are taken to be independent 
of position within the disc. 
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4.2. Data fitting 

In the process of extracting the disc geometric parameters 
from the observations, we first attempt to fit the I-, H- and 
K-band images simultaneously with a single disc model. 
This is presented in 14.2.11 We use flmax, 'Tidustdisc, P-, a 
and the inclination as free model parameters. The I band 
image suggests a maximum disc radius of 8 arcsec, which 
corresponds to 1120 AU at a distance of 140 pc. There 
must also be an inner disc radius, which is typically at a 
few tenths of an AU. These radii are poorly constrained 
by our I-, H, and K-band images; we therefore fix them at 
the given values (inner radius is set to 0.5 AU), which will 
have little effect on our results. Since we do not implement 
background illumination into our model, the appearance 
of the dark lane and the counter nebula is expected to be 
slightly different in the model. 

We computed about 12000 models covering a wide range 
of the free parameters. A handful of viable solutions remain 
after the image-fitting process (in sec. I4.2.2p . We explore 
the relative merits of these few solutions by calculating their 
respective SEDs (in sec. l4.2.3"l) . and I-band polarisation map 
(in sec. I4.2.5P for comparison with the data. 



4.2.1. Image processing and extraction of observables 

To permit a direct comparison of the observed and mod- 
elled images, observational effects must be added to the 
model images. First of all, the pixel size of the model out- 
put is therefore chosen to match the scale of the observa- 
tion (O'.'200/pixel in the optical, O'.'211/pixel in NIR). The 
model images must be convolved with the observed PSF 
and scaled to the same peak value as the observed images 
which allows an addition of Poisson and background noise 
to yield the same signal-to-noise ratio. 

Having two comparable images, it is now possible to 
quantitatively compare the model and observed images. To 
avoid potential biases induced by asymmetries in the disc 
structure and the ignorance of the background illumina- 
tion of the object and to circumvent the intrinsically noisy 
nature of Monte Carlo images, we decided to use a "qual- 
ity metric" based on geometrical observables of the neb- 
ula rather than an image-wide pixel-to-pixel goodness-of- 
fit estimate. These observables, which are detailed below 
and illustrated in Fig. [5] are chosen to describe the width 
of the nebula, its triangular shape and the counter neb- 
ula behaviour relative to the northern nebula. They are 
extracted with the same routines in the observed and mod- 
elled images, after further smoothing the images with a 
gaussian (width ct = 1 pixel). We first define the outer 
"ridges" of the triangular nebula. To allow this, we calcu- 
late contour lines at the fiux levels /„ = /max • 2^" with 
n = 1,2, ...,7V, while N defines the last contour line at 
flux level higher than the noise. The ridges (labelled with 
ai and 02 in Fig. [5]) are then defined by two fines that cross 
those points on each contour that are most distant from 
the central peak. We derive these lines as follows: We rep- 
resent the contours in polar coordinates as Vn {(t>) where r is 
the distance from the central peak and cj) is the angle from 
North. These polar functions are then normalised and aver- 
aged {r{(j>) = n-^J2n'^n{4')/ Jo'' rn{(l)')d(l)'). The direction 
of the two lines are described by the two values of (j) where 
r(0) has a maximum. Thereafter, the opening angle S be- 




Fig. 5. Definition of geometrical observables w,a,b,c and 
S. The contour lines are shown for the levels /max ■ 2~" 
n = 1 ... 10. The two ridges and the symmetry axis are 
drawn as dashed lines. Their intersection points with the 
contour lines of the flux level / — /max • 2^* are shown as 
boxes and 1=1 ■ CTnoiso as circles, respectively. The counter 
nebula peak is shown as a star. 



tween the two ridges was measured. The symmetry axis of 
the nebula corresponds to their bisecting line. 

The width of the nebula is defined as the distance w 
between the two intersection points of the ridges with the 
contour at the level of a signal-to-noise ratio equal to 1 
(shown in Fig. [5] as circles). 

The next observable describes quantitatively whether 
the nebula looks triangular or round. This observable, 
which we call the triangularity, t, is determined for each 
contour level. In Fig. [5] an example is shown for the con- 
tour level n = 8. The two ridges and the symmetry axis 
are intersected with the contour line (shown as squares in 
the figure). The distance between the intersection point of 
the symmetry axis with the contour line and the triangle 
spanned by the peak and the two intersection points of the 
contour line with the ridges is called b. The triangularity is 
then defined as 



a(l - cos((5/2)) 



(3) 



where a is the average of ai and 02. This parametrisation 
was chosen to describe a perfectly triangular nebula for 
t = and a perfect circle for t = 1. For the model-to- 
data comparison, we fitted the most triangular contour line, 
which describes best the shape of the outer nebula. 

Finally, the counter nebula is described by the peak-to- 
peak distance c, and the contrast is defined as the ratio 
between the maximum brightness of the nebula and the 
counter nebula. 



4.2.2. Image fitting 

For a given set of free model parameters (a, /3, ...), we use 
the following quality metric to estimate whether the model 
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Fig. 6. Examples of projections of f^ (eqn. |3|) to different 
model parameter planes. The contour is plotted in a loga- 
rithmic scale. White areas correspond to minimal values of 
/^. The colour range is identical in all four plots. The black 
cross shows the location of the best model, after the SED 
fitting. 

actually reproduces the three images of the object. This 
metric is a pseudo-x^, defined as follows: 



f-E 



\,i 



{ ^raodel „> 



,obs\2 
i,x) 



• 9i,X 



(4) 



where o™^^"^'^' and o°^ are the observables with index i (i.e. 
w, t, c, S, contrast) for a given photometric band A (I, H, 
K) of the model image and observation image, respectively. 
The normalisation with Wi^x, which is determined by the 
standard deviation std{{o™'^'^'^^ — of'^Y) over the whole ex- 
plored model parameter space, allows an equitable compar- 
ison of different quantities. Additionally, each observable is 
weighted with gi^\. These weights were fine-tuned by hand 
to equalise the influence of each individual parameter. This 
weighting is important to allow for a wavelength-dependent 
contrast parameter (e.g. the counter nebula is not visible 
in the H and K band, therefore the parameters c and the 
contrast are useless for these bands). Additionally the use 
of some parameters provides stronger constraints towards 
a good fitting solution than others (e.g. the nebula width 
and triangularity in comparison to the opening angle). The 
weights are listed in Table [TJ 



Table 1. Weights gi^\ of the pseudo-x^ function ((4]) 





d- 


It; 


t 


c 


Contrast 


I 

H 
K 


1/9 
1/7 
1/7 


3/9 

3/7 
3/7 


3/9 

3/7 
3/7 


1/9 




1/9 





The exploration of the parameter space was performed 
iteratively. We first ran a series of models with only 2 vari- 
able model parameters. Once the best model in the series 



Fig. 7. Comparison between observations (left) and best 
fitting model (mid, see Table [2]) in I-band (top), H-band 
(middle) and K-band (bottom). For demonstration pur- 
poses, the right column shows the images of the best fitting 
model but with doubled dust disc mass (3 ■ IQ'^^Mq). 



was identified, we fixed these model parameters and cre- 
ated a new series by varying two other parameters, and 
so on. The range of the explored parameter space is listed 
in the second column of Table [2l The first selection of the 
parameter range was chosen arbitrarily but reasonably. If 
the model behaviour suggested possible solutions outside 
the scanned range, the exploration was extended until the 
model diverged significantly from the observation. Among 
all these models, only a few provided good fits simultane- 
ously to the images of IRAS 04158-1-2805. In addition, we 
note that fitting the images at each wavelength indepen- 
dently leads to models of similar parameters and quality. 
The combination of all images in a single fit allows, how- 
ever, to narrow down the list of possible good models. 

Fig. [6] shows the map of /^ as a function of pairs of 
free model parameters: Dust disc mass vs. inclination (H^), 
flaring coefficient /3 vs. inclination ^p), dust disc mass 
vs. maximum grain size ([6];), (3 vs. scale height hp (l6]i). 
These plots demonstrate the complexity of the search for 
the minimum since some model parameters are strongly 
correlated (shown, e.g., in Fig. [6^ and b) or show numer- 
ous local minima (e.g.. Fig. [6ji). While Fig. [61; implies a 
smooth topology for the Mdust vs. amax plane, the mini- 
mum found at Mdust — 1-5 • IQ^'^Mq and Omax = 2/im is 
only valid for the inclination of i = 63°. At an inclination 
of ? = 61° we find a minimum with a similar goodness-of- 
fit with Mdust — 2.5 • 1O~*M0 and Omax — 2.5^m. This 
implies that we cannot find a unique best model by fit- 
ting images only. The synthetic images for our best-fitting 
model, found after the SED fitting process, is presented in 
Fig. [7] alongside the observed images. 

We estimated the valid model parameter range by us- 
ing one good-fitting parameter set and changed only one 
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parameter until the /^ function showed a significant in- 
crease or until an eyeball comparison could clearly detect 
an unsatisfying model. This method does not include an 
estimation of the general solution range but gives an im- 
pression of the validity of the free model parameters at a 
certain minimum of /^. 

Additionally we had to set an artificial upper limit for 
the valid range of a,nax: Fig- Efc shows a valley for Mdust = 
1.5 • 1O~^M0 and Cmax > 2/im. This is a consequence of 
comparing the modelled with the observed images only up 
to the K-band. Therefore, the fit is only mildly sensitive 
to larger grains, and we set amax = 2/i?7i with no further 
consequences for the solution. 

Fitting the images for each wavelength individually de- 
livers only very small variations in Omax and disc mass and 
is consistent with the multi-wavelength solution within the 
valid model parameter range. 

4.2.3. SED fitting 

To refine the model selection we fit the SED. The mod- 
elled SED includes photospheric emission, scattered light 
and thermal emission by assuming passive heating of the 
dust in the disc (Pinte et al. 2006). We use a 3000K spec- 
trum from Baraffe et al. (1998) in agreement with a spectral 
type of a low mass star. For each of the best image-fitting 
models there are only 2 extra free parameters to the SED 
fitting: the value of foreground extinction Ay and the bolo- 
metric luminosity of the central source. The latter is a free 
parameter since we do not see the source directly. The fore- 
ground extinction Ay describes the material between the 
observer and the outer limit of the star-disc system and is 
added to the model after the radiation transfer calculation. 

It turns out that only one model is a simultaneous good 
fit to the images and to the SED. We consider it our best 
model in the following (see third column of Table [2]) . The 
SED for this model is presented in figure [31 It includes 
Ay = 0.5 mag and a stellar luminosity of ^^0.4 Lq. 

Photometric variations are observed for 
IRAS 04158-1-2805. It results in a large spread in the 
SED data at a given wavelength in the optical and NIR 
(e.g., ~1.5mag in the I-band between Strom & Strom 1994 
and Luhman 2000). The calculated SED agrees well with 
the data points and falls within the range of observed 
variations. The estimated extinction is lower than the value 
estimated by Strom & Strom (1994; Aj=0.7), Luhman & 
Rieke (1998; Aj=1.5), and Luhman (2000, Aij = 1.5) in the 
NIR. However these estimations did not take into account 
the fact that the object's colours are strongly modified by 
the scattering in the circumstellar material. We therefore 
believe that our estimate of the foreground extinction is 
more robust than previous estimates. 

4.2.4. Acceptable parameter range 

The parameters of the best-fitting model are described in 
Table [5] under " best model" . To define the " acceptable 
range" (see fourth column of Table ^ we calculated series 
of images from the best model in which only one parameter 
at a time was varied from its best fit value. This allowed us 
to define a possible range of solutions, i.e., solutions with 
values of p only moderately higher than the best model 
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Fig. 8. Difference of observed and modelled polarisation 
percentage (top) and polarisation direction {hottom) along 
the ridges (black) and the symmetry line (grey). 



matches with the observations by an eyeball comparison. 
These ranges give a feeling of how tightly each parameter 
is constrained. 



Table 2. Disc model parameters for the best-fitting model 
and their acceptable range. 



Parameter 


Explored 


Best 


Acceptable 




range 


model 


range 


tSmax \^J.m\ 


0.5 - 4.0 


2.0 


1.6 - 2.8 


Mdustdisc [IO-^Mq] 


0.25 - 3.0 


1.5 


1.0 - 1.75 


i?out[^C/] 


not varied 


1120 


- 


RiMU] 


not varied 


0.5 


- 


P 


1.0- 1.3 


1.20 


1.15- 1.21 


a 


-2.0- -0.1 


-1.0 


-1.5- -0.5 


ho{r ^ 50AU)[AU] 


5.0 - 12.0 


8.0 


7.8 - 8.2 


ind.[°] 


0.0 - 90.0 


62.7 


62.1 - 63.25 



(p/p 



< 



125%) or images that showed no obvious mis- 



4.2.5. Using polarisation to confirm the model 

For the best-fitting model we compared the polarisation 
map with the data. Observational effects are applied to the 
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model polarisation output in a similar manner as for the 
direct images. 

To compare the model with the data, the polarisation 
values are extracted along the ridges and the symmetry axis 
of the Northern nebula (see Fig. 2] left and middle panel) 
and compared as a function of distance from the central 
source (right panel) . The difference of the observed and the 
modelled polarisation rate and orientation can be found in 
Fig. [8l The observed polarisation rates are reproduced rea- 
sonably well along the ridges by the model {AP/P w 0.3) 
except on the central peak where the model predicts only a 
1.2% polarisation instead of the observed 3.3%. The large 
difference in polarisation orientation close to the central 
peak is an artifact from the exact peak placement and has 
no physical meaning. The upper limit of the observed polar- 
isation is also well reproduced (^30%). Along the symmetry 
axis of the nebula the general trend is also correct, rising 
from central source to edge, but the model polarisation ap- 
pears to be systematically lower by a few percent compared 
to the data. Nonetheless, the trends are well reproduced as 
well as the observed maximum polarisation rates. 

Table 3. Extracted observables of the observation and the 
best-fit model 



Observable 


Band 


Observation 


Model 




I 


15.8 


14.9 


w [arcsec] 


H 


14.1 


14.0 




K 


10.8 


12.1 




I 


129 


136 


sn 


H 


129 


133 




K 


136 


139 




I 


0.31 


0.33 


t 


H 


0.44 


0.33 




K 


0.53 


0.23 




I 


340 


300 


Contrast 


H 


>1100 


- 




K 


>8700 


- 


c [arcsec] 


I 


4.76 


4.22 




0.7fim 


0.18" 


0.25 




0.9fim 


0.90' 


1.11 




1.25/im 


1.96^-2.09=* 


2.26 




1.6pim 


2.24'',2.9^3.423 


2.85 


Flux 


2.22^m 


3.06'',3.46^4.3l3 


3.61 


[lO-^VKm-^Hz-i] 


3.8^m 


1.19^ 


4.75 




12pm 


6.50=* 


7.19 




25nm 


8.63^ 


6.93 




60fim 


7.39^ 


6.63 




ISOOfim 


0.025^ 


0.008 


min. Pol. 


I 


3.3% 


1.2% 


max. Pol. 


I 


30.9% 


30.3% 



References: 1 - Strom & Strom (1994); 2 - Luhman & Rieke 
(1998); 3 - Kenyon et al. (1990); 4 - Luhman (2000); 5 - 
Motte & Andre (2001); 6 - This paper, see Sect.|321 



5. Discussion 

5.1. Reliability of the models 

We were able to fit a dust disc model to the observations 
of IRAS 04158-f 2805 that matches the I-, H- and K-band 
images as well as the SED from visual to far-infrared and 
the I-band polarisation rate at different object positions. 



Since our model uses simple assumptions, such as homoge- 
neous spherical grains, power-law distributions (grain size, 
surface density, flaring) and no dust settling, these results 
provide a useful insight on the geometry of the circumstellar 
environment. 

We fitted the observation well in terms of the width of 
the nebula, opening angle and peak-to-peak distance and 
brightness contrast (see Table [3] for all observables). The 
triangularity t is well fitted for the I- and H-band but for 
the K-band the model is too triangular. At the centre of 
the object, the model predicts a polarisation rate which is 
slightly below the observed value. While this probably in- 
dicates that our grain model needs refinement (either in its 
size distribution, composition or grain shape), the former 
may indicate that we need a more complex geometry either 
for the disc or for the emitting source. For instance, K- 
band emission from the inner parts of the disc might be an 
important source of photons, which is neglected in our cal- 
culations. To check this possibility we calculated a K-band 
model image that includes the disc emission as a source for 
scattered photons for our best-fitting model. It turns out 
that the star is so cool that, at the 0.5 AU assumed inner 
radius, the dust reaches a temperature of only 400 K. As 
a consequence, roughly 99% of all K-band photons emitted 
by the system come from the star and, therefore, the re- 
sulting images are unchanged if we include or neglect the 
disc emission. Vertical settling could also play a role, as 
well as the presence of a remnant halo that could be re- 
sponsible for the roundish 4-5" structure best highlighted 
at K-band, where the extinction of the material located in 
front becomes much lower. 

By comparing the I-band observed and modelled im- 
ages, the shape of the counter nebula and the dark lane 
between the two bright areas do not coincide exactly. The 
model did not take into account the absorption of back- 
ground light in the dust lane. The artificial noise of the 
model images was produced by adding Poisson noise to 
the scattered light images and by superposing position- 
independent Gaussian-noise. This does not reflect the real 
nature of the background since IRAS 04158-1-2805 is il- 
luminated by the large reflection nebulosity in the back. 
Therefore, we may expect a difference of the observation 
and the models at those wavelengths where the background 
light dominates the noise. In the H- and K-bands this dom- 
inance is not visible in the observed images. 

The model underpredicts the millimeter flux by a fac- 
tor of 3 but uncertainties are large regarding dust opac- 
ities in this regime, typically by a factor of 5. Also, the 
mm-data were obtained with an 11 arcsec beam and may 
suffer from background contamination, hence overestimate 
the true flux. Higher resolution millimeter data is needed 
to investigate this discrepancy further. 

Nevertheless, one single model can describe accurately 
the scattered light images in three wavelength bands, the 
SED and the polarisation map of IRAS 04158-1-2805. The 
method applied here is therefore generally valid and promis- 
ing to find parameters for relatively simply structured dust 
discs. 

5.2. Gas to dust ratio 

From our model, we can easily derive the column den- 
sity adust by integrating the dust density structure along 
our line of sight to the central source. We find (Tdust = 



10 
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From Nh we obtain a gas column 
X 10^^ g cm^^ by assuming inter- 



i2t\% X 10-"* g cm-2. 

density Ggas = T-'^ttl 
stellar abundances (Morrison & McCammon 1983). This 
provides an estimate of the total gas-to-dust ratio in a pro- 
toplanetary disc, along the line of sight that grazes the disc 
top layers in this case. To our knowledge this is the first 
time it is obtained directly using Nh rather than Nco'- 



gas _ r,r,p, + 170 



<^dust 



(5) 



While uncertainties on this ratio are still substantial, 
this illustrates a robust method to derive dust-to-gas ra- 
tios in protoplanetary discs. Although the value we find is 
compatible with the standard value of 100 usually assumed 
in this context, our result might suggest a slightly higher 
value in the top layers of the disc. This is a potentially im- 
portant result that calls for further investigations in more 
sources. 

In any case, we do find an upper limit to the gas-to-dust 
ratio (albeit not a lower limit) of « 560(2ct), which provides 
clear evidence that the disc cannot be strongly (e.g. by an 
order of magnitude) depleted of dust along the line of sight 
to the star. Such dust depletion would be expected if strong 
settling of dust toward the disc midplane had occurred. 

5.3. IRAS 04158+2805: a classical T Tauri star with a large 
disc? 

We have shown that we can well reproduce all observational 
properties of IRAS04158 with a simple model of a nearly 
edge-on disc, i.e., that of a Class II source without sub- 
stantial circumstellar envelope. This result confirms prior 
classifications by e.g. Park & Kenyon (2002) or Kenyon 
& Hartmann (1995). The flat SED of IRAS 04158-^2805 
can be interpreted in terms of a high inclination to the 
line of sight instead of invoking a highly embedded source, 
as a protostar would be. This work demonstrates the im- 
portance of understanding the circumstellar geometry to 
assess the nature of the central source. Here, an edge-on 
disc blocks our direct view towards the central object, flat- 
tening the shape of the SED. It is known that a classical 
T Tauri star can have a rising near- to mid- infrared SED. 
But only a few mimic an embedded source, as is the case for 
IRAS 04158-^2805. Most sources such as HH30 (see, e.g.. 
Wood et al. 2006), HK Tau B, or HV Tau C, do have a 
declining near- to mid-IR SED, and the second peak from 
thermal emission is found only at longer wavelengths. 

Interestingly, IRAS 04158+2804 is just above the sub- 
stellar limit, based on its spectral type. If our model is 
correct, it hosts a large massive disc implying that some 
of the lowest-mass T Tauri stars, at least in Taurus, can 
be surrounded by ^lOOOAU-radius discs. This object likely 
formed from the collapse of a small pre-stellar core and 
is very unlikely to have undergone any violent dynamical 
interaction, such as an ejection from an unstable multiple 
system (e.g. Reipurth & Clarke 2001). This result adds sup- 
port to the idea that the ejection scenario is not the only 
mode to form VLM T Tauri stars and brown dwarfs. 

Also, considering the probable mass of the central 
star r^ 0.1 — O.2M0, the total dust-disc mass we derive 
(1 — 2-10~'*Mo) and the gas-to-dust ratio we derived (^^220, 
assuming it is true throughout the disc and not only along 
our line of sight), this implies a total disc/star mass ra- 
tio of ~0. 1-0.2. In other words, the disc is close to the 



limit against gravitational instability. Because collapse in 
these gravitationally unstable discs is one suggested mode 
for planet formation, it is of great importance to study 
IRAS 04158+2805 further: The estimation of both the gas 
and the dust masses need to be refined. The gas mass esti- 
mation can be improved by much deeper X-ray observations 
while the dust mass can be obtained better with models 
witch use more free parameters and which are directly fit- 
ted on the images, not just a handful of observables. 

6. Conclusion 

In this paper we presented a multiwavelength study of 
IRAS 04158+2805 and its circumstellar environment. We 
modelled the shape and brightness profiles of the reflection 
nebulosity in three optical and NIR bands (I-, H- and K- 
band) with MCFOST, a Monte Carlo polarised radiative 
transfer code, and found a good agreement between the 
model and the data. Many parameters of the final model are 
well constrained (e.g., inclination, Omax, -^dustdisc) while a 
few remain poorly determined. The scale height h and (3 
are degenerate but the pair h-/3 is relatively constrained. 

The disc model parameters used to match the data are 
listed in Tabled They are 2fim for the maximum dust grain 
size, 1.5-lO^^Af0 for the dust disc mass, 1120 AU for the 
outer disc radius, 1.2 for the exponent f3 of the flaring law, 
-1.0 for the exponent a of the radial dust density law, a 
scale height of 8AU at a radius of 50 AU and an inclination 
of 62.7°. 

The best model fits the observed SED well. However, it 
falls slightly short of producing the right amount of 1.3mm 
continuum flux by about a factor of 3. The model also re- 
produces reasonably the observed I-band polarisation be- 
haviour along the symmetry axis and ridges of the Northern 
nebula. 

Combining dust disc models with X-ray spectroscopy 
allowed probing the gas-to-dust ratio in the disc. We found 
a value of 2201|^5q which is compatible with the ISM value 
and the value generally assumed in protoplanetary discs. 

According to its spectral type, IRAS 04158+2805 has a 
mass slightly above the substellar limit. Clearly, stars with 
such a low mass keep being formed by accretion from a 
circumstellar disc whose properties do not seem to differ 
significantly from those of their more massive young coun- 
terparts. It would be interesting to push the search further 
for disc images around less massive objects, located well 
into the substellar regime. 
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